function section = build_side_section(grids,ddc_grids,f)
% side_v_x = build_side_section(grids,ddc_grids,f)
%
% Build a side section corresponding to the implicit surface
%  f(x) = 0
% The section is oriented from negative values of f toward positive
% values.
%
% Note: this function requires the data structures produced by the
% mod_grid module. So it is necessary to run some fortran the code once
% before using this subroutine. The number of partitions used in the
% fortran run, however, does not matter.
%
% grids : cell array with all the grids to be sampled
% f     : function pointer

% Algorithm:
% 1) identify all the elements with barycenter such that f(xb)=<0
% 2) loop over the sides of such elements: all the sides connected to an
%   element with f(xb)>0 are included in the section, with the
%   appropriate orientation.
% 3) for boundary sides, the barycenter of the side is used to check
%   the sign of f.

 side_v_x = [];
 for i=1:length(grids)

   grid     = grids{i};
   ddc_grid = ddc_grids{i};

   el_pos = zeros(1,grid.ne);
   for ie=1:grid.ne
     el_pos(ie) = (f(grid.e(ie).xb)>=0);
   end

   ns = 0;
   lside_v_x = zeros(grid.m,grid.d,grid.ns);

   for ie=1:grid.ne
     if(not(el_pos(ie))) % negative
       for is=1:grid.d+1
         skip_side = 0;
         ie2 = grid.e(ie).ie(is);
         if(ie2>0) % internal side
           pos2 = el_pos(ie2);
         elseif(ie2==-ddc_grid.ddc_marker) % ddc side
           id = ddc_grid.ns(2,ddc_grid.gs(2,grid.e(ie).is(is)));
           in = ddc_grid.ns(3,ddc_grid.gs(2,grid.e(ie).is(is)));
           ie2 = grids{id+1}.s(in).ie(1);
           x2 = grids{id+1}.e(ie2).xb;
           pos2 = (f(x2)>=0);
           % Periodicity makes things more complicated here. To reduce
           % the problem, we check that the two elements lay on the
           % opposite sides of the side.
           nn = grid.e(ie).n(:,is)'; % normal
           n1 = nn * ( grid.e(ie).xb - grid.s( grid.e(ie).is(is) ).xb );
           n2 = nn * (      x2       - grid.s( grid.e(ie).is(is) ).xb );
           skip_side = (n1*n2)>0;
         else % boundary side
           x2 = grid.s( grid.e(ie).is(is) ).xb;
           pos2 = (f(x2)>=0);
         end
         if(and(pos2,not(skip_side))) % found a section side
           ns = ns + 1;
           for j=1:grid.d
             lside_v_x(:,j,ns) = grid.v( ...
                 grid.s( ...
                   grid.e(ie).is(is) ...
                 ).iv( ...
                   grid.me.pi_tab( grid.e(ie).ip(is) ).i(j) ...
                 ) ...
               ).x;
           end
         end
       end
     else % positive
       % These sides can still be included if they are boundary sides,
       % since in such acase there is no risk of duplication
       for is=1:grid.d+1
         ie2 = grid.e(ie).ie(is);
         if(ie2<0) % boundary side
           if(ie2!=-ddc_grid.ddc_marker) % not ddc
             x2 = grid.s( grid.e(ie).is(is) ).xb;
             if(f(x2)<=0)
               ns = ns + 1;
               for j=1:grid.d
                 lside_v_x(:,j,ns) = grid.v( ...
                     grid.s( ...
                       grid.e(ie).is(is) ...
                     ).iv( ...
                       grid.me.pi_tab( grid.e(ie).ip(is) ).i(j) ...
                     ) ...
                   ).x;
               end
             end
           end
         end
       end
     end
   end

   side_v_x = cat(3,side_v_x,lside_v_x(:,:,1:ns));
 
 end

 section = struct( ...
    'section_type','side_section', ...
    'section_name','<your section name>', ...
        'side_v_x',side_v_x ... 
 );

return

